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Charged systems interacting via Coulomb forces can be efficiently simulated by introducing a local, diffusing 
degree of freedom for the electric field. This paper formulates the continuum electrodynamic equations corre- 
sponding to the algorithm and studies the spectrum of fluctuations when these equations are coupled to mobile 
charges. I compare the calculations with simulations of a charged lattice gas, and study the dynamics of charge 
and density fluctuations. The algorithm can be understood as a realization of a mechanical model of the ether. 



I. INTRODUCTION 

Molecular dynamic simulation of charged condensed mat- 
ter systems is slow and difficult 1 . In standard methods, such 
as optimized Ewald summation 2,3 , fast multiple methods 4 - 5 
or the fast Fourier transform 6 , extensive and time consuming 
bookkeeping is needed because of the range of the Coulomb 
interaction. This bookkeeping often scales badly when imple- 
mented on modern multiprocessor machines which are used 
in the simulation of the largest systems. Naive Monte-Carlo 
methods are particularly inefficient since the motion of a sin- 
gle particle in an N particle simulation requires the recalcula- 
tion of the Coulomb interaction with all other particles, lead- 
ing to a complexity in O(N) for an update in the position of a 
single particle 

Recently, 7 a local algorithm with complexity scaling as 
O(l) per update was introduced for the Monte-Carlo simu- 
lation of charged particles. In this algorithm an auxiliary elec- 
tric field E is coupled to the charge density. The dynamics 
of E are chosen so that the equilibrium distribution is deter- 
mined by the Coulomb interaction. Due to the locality of the 
algorithm the method is trivial to implement on multiproces- 
sor machines. In this paper we study the dynamics of the al- 
gorithm in order to understand the relaxation processes and 
time scales involved in a typical simulation. Simulations are 
performed on a model of a charged lattice gas to demonstrate 
the diffusive propagation of charge and density fluctuations. 

The algorithm is based on implementing Gauss's law 

div E = p/e Q (1) 
in the equivalent integral form 

J-E.dS = q/e (2) 

as an exact dynamic constraint on the Monte-Carlo algorithm. 
Here p is the charge density and e the dielectric constant and 
q the charge enclosed by the surface of integration in eq. (|2jl. 

The system is discretized by placing charged particles on 
the vertices of a cubic lattice, {i}. The electric field E^j is 
associated with the links of the lattice. The simulation 

starts with Gauss's law satisfied as an initial condition. There 
are two possible Monte-Carlo moves: Firstly, fig. Q, we dis- 
place a charge, e, situated on the leftmost lattice site, 1, to the 
rightmost site, 2. The discretized constraint is 

^2 E i,j = e i/ e Q ( 3 ) 
i 
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FIG. 1 : Motion of a charge e from 1 to 2 is associated with a modi- 
fication of the field on the connecting link: Ei,2 — > £1,2 — e/eo. 
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FIG. 2: Modification of the angle Q Z A leads to modified values of 
the field between on the links {1, 2}, {3, 2}, {4, 3} , {4, 1}. The 
field in the y direction associated with the link is given by the 
difference of the angles associated with the two plaquettes, Q\ and 
Qg, which are aligned in the x direction. 



with ei the charge at the site i. The sum in eq. corresponds 
to the total electric flux leaving the site i. The constraint is 
again satisfied if the field associated with the connecting link 
is updated according to the rule 2 — * Ei^ — e/eo. Secondly 
we update the field configurations, fig. (0,, by modifying the 
four field values of a single plaquette by a pure rotation; Ei^ 
and E^i increase by an increment A whereas E4.3 and £3.2 
decrease by A so that at each vertex the sum of the entering 
and leaving fields is unchanged. The basic dynamic degree of 
freedom in the second update is a circulation or rotation, 0, 
associated with each plaquette of the network. 

In the first section of the paper we formulate the continuum 
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limit of the evolution equations and show that they lead to 
diffusive evolution of the electric field. We then couple the 
electric field to a mobile gas of charged particles and compare 
the solutions of the coupled plasma equations to simulations. 
Finally we show that the equations are closely related to the 
Maxwell electromagnetic theory. 



II. DIFFUSIVE ELECTRODYNAMICS 

A. Fundamental equations 

We start with a simple example to motivate our derivation 
of the effective large scale equations obeyed by the electric 
field: a single particle diffusing in a harmonic potential with 
energy U = Kx 2 /2. The equation of motion is found by 
taking the derivative of the energy with respect to the dynamic 
coordinate x and then balancing the resulting force against a 
relaxation process linear in the velocity 



(4) 



where f(t) is an external forcing term. £ the inverse mo- 
bility sets the characteristic time scale of the relaxation of 
x and is a function of the step size of the Monte-Carlo trial 
moves. A first order (in time) algorithm, such as Monte- 
Carlo, for the simulation of a particle in such a potential (per- 
formed near equilibrium with small step sizes) is essentially a 
discretized realization of this stochastic differential equation. 
The Langevin description is completed by specifying /, so as 
to obey the fluctuation dissipation theorem. 

We now turn to the equations for the electric field. Firstly 
we examine the field in absence of current then we shall find 
the coupling of the field to external sources. The discretized 
energy of the electric field is given by 
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(5) 



The basic variables in the dynamics of the field are, however, 
not E but rather the rotational degrees of freedom which are 
updated independently at each time step. The conjugate force 
acting on this variable is a torque. We define on each plaquette 
in the {a;, y} plane the angular variable Z . We group the an- 
gle corresponding to the three possible plaquette orientations 
into a vector 0. The physical field E is sensitive to deriva- 
tives of 0: Studying fig. (0 one sees that E\ t 2 is given by the 
x variation of the z component of 0. This we recognize as 
part of the curl operator acting on the field 0. The complete 
relation between a local variation in E and a local variation in 
is thus 



<5E = curl 6® 



(6) 



where in the discretized model derivatives are to be under- 
stood as finite differences with E living on links and on 
plaquettes. 



Modification of the circulation Q Z A in fig. (|2ji by A gives 
rise to a new contribution to the energy of the plaquette. 



U = J ( ( ^4,1 + A)' + {Ei,2 + A) 2 4 
( £4,3 - A) 2 + (£3,2 - A) 2 ) 



(7) 



Taking the derivative with respect to A we find a torque, C, 
acting on the circulatory degree of freedom which is a dis- 
cretized version of the curl of the field. 



_ dU 

Cz ~~dA 



e ( (£4,1 - £3,2) + (£1,2 - £4,3)) (8) 

(9) 
(10) 
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Again the three components of the torque live on the plaque- 
ttes together with the angle 0. By analogy with eq. @ the 
evolution equation for the angle © is 
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-cq curl E 



(11) 



linking a velocity to the conjugate force. Eventually a stochas- 
tic force should also be added into this equation but we shall 
not need it in what follows. 

Consider now the evolution of the field in the presence of a 
current. Take a network in which charges are present at every 
vertex and displace every charge to the right as in fig. {TJ. 
Then every bond in the x direction is modified by — e/eo 
where e is the charge displaced even though the local charge 
density is unchanged. If we displace the charges at a constant 
rate we have the evolution of the field due to the source as 



9E S 



dt 

Combining eqs. d6ll 2t we find 
<9E 



-J/e 



at = J/e ° + curl ~dT 



(12) 



(13) 



This equation is clearly analogous to the Maxwell equation, 
<9E 



co- 



dt 



curl H 



(14) 



if we write that H = €07^. We shall see later that H is indeed 



closely related to the magnetic degrees of freedom of classical 
electrodynamics. From eqs. (1 1 01 1 11131 



<9E 

— = (e V 2 E - grad p) /£ - J/e 



(15) 



where we have used the standard identity curl curl E = 
( grad div E — V 2 E) and Gauss's law. Equation J15I is the 
main result of this section giving a diffusive propagation law 
of the electric field in absence of external charges and cur- 
rents. 
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In the static limit both the current and the time derivative 
of eq. d!5t vanish. We find the same equation for the electric 
field 



V 2 E = grad 



(16) 



as is found by applying the operator (— grad ) to the Pois- 
son equation in conventional electrostatics. When we take the 
divergence of eq. ( 1151 we discover that 



dt £ 



(divE-p/eo) = 



(17) 



Again we see that Gauss's law is implemented in the method 
as an initial condition. Note that for this to be true we require 
that the Langevin noise associated with eq. (1151 does not in 
itself destroy the conservation law. It is thus the curl of some 
vector field 8 . 



B. Relation to Potentials 

In our earlier paper 7 we showed that the electric field could 
be calculated from a scalar potential and a vector potential 
Q with the relationship 



E 



grad i 



curl Q 



In Fourier space we can write this equation as 



E(k) 



-ik</> + ik A Q 



(18) 



(19) 



The second term of this expression is perpendicular to k so 
that there are two transverse degrees of freedom in the Q field, 
corresponding to two independent polarization states; the lon- 
gitudinal component of Q is projected out and does not con- 
tribute to the electric field. We can consider that the field is 
due to a static longitudinal potential plus transverse photons. 

If we take the time derivative of eq. i ll 8b we can com- 
pare with eq. Jl 3b . The term curl is purely transverse, 
whereas —J contains both longitudinal and transverse com- 
ponents. Thus we conclude that 



and 



curl Q = curl — J t /en 



grad <j> = Ji/e 



(20) 



(21) 



where J t and J/ are the transverse and longitudinal compo- 
nents of the current. In general these are non-local relation- 
ships since the projection involves a passage via Fourier com- 
ponents. Such non-local relationships between potential and 
field are normal in the Coulomb gauge 9 . 



C. Phenomenological Dynamics of a Two Component Plasma 

In this section we couple the diffusive evolution equation 
for the electric field eq. (1151 to the equations describing a two 



component plasma and study the relaxation phenomena and 
time scales that are to be expected when using the algorithm 
to simulate dense charged systems. 

The equations of conservation and linear response give the 
following equations for the charge degrees of freedom: 



dc+ 
~df 

~df 
J+ 
J 



= -divJT 



div J" 



D grad c 1 
-D grad c" 



e^c + E 
e/ic~E 



(22) 



Where c + and c~ are the number densities of of the positive 
and negative charges, ±e. J ± are the number current den- 
sity. From these equations we find the equations obeyed by 
the total density c and the charge density p. We note that the 
diffusion coefficient D and the mobility, p are linked by the 
Einstein relation D = ksT p. 

Taking the sum and difference of the equations ( 1221 we find 



dc 

dt 
dp 

dt 

J 



DX7 2 c 



div J 



= eVcoE — D grad p — e 2 cop grad (j> e {t) (23) 



where we have linearized the equations about a mean density 
Cq. J = eJ + — eJ is the electric current density. <f> e (i) is 
an externally imposed potential that we shall use to calculate 
charge-charge correlation functions. Substituting eq. ( 1231 for 
J in ( 1151 we find 



<9E 

~dt 



= e /£V 2 E- e 2 M c /e E (24) 
+ (D/e - 1/0 grad p + e 2 c /x/e grad cj) e 



We analyze this equation by treating separately the longitudi- 
nal and transverse fluctuations. Take the curl of eq. ( 1241 to 
find that the transverse components of E decouple from the 
charge density. In Fourier space we find the dispersion law 



(ioj + e /£ q 2 + e 2 /ic /e ) q A E = 



(25) 



In the absence of charges the mode is diffusive but the pres- 
ence of a finite charge density gives a gap in the spectrum. 

Consider now the equations for the field eq. ( 1241 : With the 
help of Gauss's law one replaces the divergence of the field by 
the charge to find 



— - DV 2 + e 2 /ic /e ) p = c Q e 2 pV 2 (j) e (26) 



which also applies to the longitudinal mode of the electric 
field 

(iuj + Dq 2 + e Vo/e ) q • E = -e 2 c pq 2 (t)e/eo (27) 
Again the spectrum has a gap as q — > 0. 
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FIG. 3: Fit of a density-density, top, and a charge-charge, bottom, 
correlation functions eq. <29> to a single exponential. 5000 particles, 
network of 25 x 25 x 25, mode q = 2-k x (2, 2, 0). Arbitrary units 



III. NUMERICAL RESULTS 
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FIG. 4: Characteristic time extracted from S(q, t) as a function of q 2 . 
Bottom curve : Density-density correlations: the mode is diffusive, 
eq, <23> . Top curve: The charge-charge correlation function has a gap 
in agreement with eq. 1261 . Selected modes between (q/2n) — 1, 
q = 2tt x (1, 0, 0) and (g/2vr) 2 = 18, q = 2tt x (4, 1, 1). 1000 
Particles, 18 x 18 x 18 mesh. Vertical axis in arbitrary units. 



A. Dynamics 

We performed simulations of a charged lattice gas to study 
the dynamics of the density and charge fluctuations. Equal 
numbers of positively and negatively charged particles with 
e = ±1 were placed on the vertices of a network which 
was simulated by the algorithm in a uniform dielectric back- 
ground. During the simulations we measured the Fourier 
transform of the particle distributions 

s(q, t) = — L e * exp(ir* (i) • cy) (28) 

i 

where the weight ej is the charge for the charge correlation 
function and is unity for the density correlation function. We 
use this information to construct the dynamic structure factor 

S(q,t) = <a(q 3 i)a(-q,0)) (29) 

The result is fitted with an exponential and the decay rate plot- 
ted as a function of q 2 in figs. (3,4). The density-density 
correlation function displays simple diffusive behavior. The 
charge-charge correlation function is characterized by a gap 
at q = 0. 

What do these dispersion relations imply for the equili- 
bration of a system of charged particles? The mass degree 
of freedom is diffusive so that a simulation equilibrates in a 
time which scales quadratically with the linear dimensions of 
the system. The charge degree of freedom is associated with 
a Green function which is also diffusive. However, the to- 
tal weight decays exponentially in time. The signal due to a 
charge fluctuation is very weak beyond the Debye length. 

Note that the parameters used in the derivation of the 
plasma dynamics are already at a coarse grained level of de- 
scription. We expect that the bare parameters are renormal- 
ized by non-linear interactions: While the eq. d!5l > is in some 



sense fundamental, containing within it the exact statement 
of Coulomb's law and Boltzmann statistics, the eqns |23i are 
purely phenomenological. An example is the mobility of a 
particle /i which in the above theoretical presentation appears 
independent of the field parameters eo and £. However con- 
sider the case of a charged particle pulled by an external non- 
electric force in the presence of a electric field which relaxes 
very slowly. As the particle moves it leaves behind it a "string" 
of electric field due to the dynamics of fig. Q. This creates a 
back force on the particle which reduces its mobility. Monte- 
Carlo moves on the field spread this string over many lattice 
sites increasing the mobility of a charged particle. Thus the 
mobility of the charged particles increases when the field re- 
laxes more rapidly. 

This effect is an explanation of the curves of fig. @ Despite 
the predictions of eqs. J23I26I > the slope of the charge-charge 
and the density-density curves are slightly different; the ef- 
fective diffusion coefficient of the charge fluctuations is lower 
than that of the density fluctuations. Slow relaxation of the 
electric degrees of freedom should hinder the motion of a sin- 
gle charged particle more than a strongly coupled, neutral pair 
moving in the same direction. 



B. Screening 

From the Poisson-Boltzmann equation it is known that 
charged systems screen. We derive this result from our dy- 
namic equations as follows: Consider eq. ( I26i for the charge 
density in the presence of a static external potential e (q). 

— q 2 Coe 2 (j) e 

= 2 , 2 — 7 — u rr , rr (30) 
q z + e z c /e Q k B ± k B ± 
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of the charge and density fluctuations, we find errors of only 
0(e p / L 1 / 2 ) per site. The high statistics curves of this paper 
were generated by using runs of length up to 5000 times the 
equilibration time. Even here the errors remained acceptably 
small. 
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FIG. 5: Plot of 1/5(5) as a function of (2n/q) 2 . The plot is linear 
as implied by eq. J 3 It . Selected modes between 2-k x (1, 1,0) out to 
2tv x (5, 5, 5). The plasma screens interactions exponentially. 5000 
charges on a network of 25 x 25 x 25 



from linear response theory the structure factor with the nor- 
malization of eq. d29i is given by 



S(q) 



1 2 
e z q z 



(31) 



where the inverse Debye length, k, is given by the standard 
expression k 2 — e 2 co/eo ksT . This prediction is checked in 
our code by plotting 1 / S(q) as a function of 1 /q 2 , fig. l|5}- 

Fig. Q should be taken as very strong evidence that our 
algorithm is behaving correctly. It reproduces one of the most 
striking features of charged systems, the exponential decay of 
the charge-charge correlation function due to Debye screen- 
ing. 



C. Numerical Stability 

In the simulations that we performed to study the dynamic 
and screening properties of the algorithm we were agreeably 
surprised by the numerical stability of the algorithm: At each 
update one makes an error e p comparable to the round off 
error of the computer. Over many time steps this accumulates 
so that Gauss's law is violated. We feared that this local error 
would rapidly become important. 

The slow propagation (in time) of numerical errors can be 
understood by consideration of eq. illl . Local fluctuations 
in the constraint div E — p/eo = spread out via a diffu- 
sion process. Since both positive and negative errors are made 
during a simulation there is a large degree of cancellation 
occurring. After a single Monte-Carlo sweep of the system 
Gauss's law is violated by 0(e p ) at each lattice site. However 
averaged over a sample with L 3 sites the average error per 
site is 0(e p /L 3/>2 ). When simulated for 0(L 2 ) sweeps the 
system comes to equilibrium under the diffusive propagation 



IV. PROPAGATIVE FIELD EQUATIONS 
A. Maxwell's equations 

In eq. d!5t we gave the equations of motion for the electric 
field obeyed in the continuum limit of a Monte-Carlo simu- 
lation. In this section we shall see how local imposition of 
Gauss's law can be used to find a propagative dynamics for 
the evolution of the electric field. We continue to describe the 
basic dynamic degree of freedom as an angular variable, 0, 
which is linked to the electric field by eq. Jl 3i . This variable 
is associated with a angular velocity, vg and moment of inertia 
Ig. The resulting second order equations will display propa- 
gating, wave like features rather than the diffusive propagation 
characteristic of eq. Jl 51 . 

Each link field E^j is the result of the rotation of a variable 
defined on the faces of the cube rotating at angular velocity 
vg . As above the torque on the rotational degree of freedom 
of a plaquette is given by C = — eo curl E. Using eq. ( 161 121 
we find the following equations obeyed by the fields: 

dv e 



Ig— — = -e curlE 
at 

<9E 

~dt 

div E 



curl Vg 

P/eo 



J/eo 



(32) 



where the differential operators are to be interpreted as the ap- 
propriate difference when acting on the lattice variables. The 
equations (1321 are a rescaled version of Maxwell's equations 
with vg playing the role of the magnetic field H. 

In order to find the coupling between particles and the vari- 
ables we are obliged to use the formalism of Lagrangian 
dynamics. Naive arguments based on energy considerations 
are ambiguous and can easily lead to wrong results. 



B. Lagrangian Treatment of Dynamics 

We shall now show how to derive the full coupled equations 
between particle motion and field. Firstly, however, we shall 
look at a simple illustrative example in order bring out the 
main formal features of constrained Lagrangian dynamics. 

Consider two gears described by the rotation angles <p and 
ip. We take these gears to have unit inertia and impose on them 
a potential energy g((p) andh(ip). The gears are in contact and 
are thus submitted to the rolling constraint 

(p + ip = (33) 

We find that the Lagrangian describing this system is simply 

£_ Jj?_ 
2 2 



£ 



g-h + A(<p + 



(34) 



6 



where the Lagrange multiplier A imposes the constraint. Note 
that we are not using the standard method of D'Alembert 
of imposing nonholonomic constraints but rather the "vako- 
nomic" method 10 in which the field A is itself considered an 
independent dynamic variable. Such methods are widely used 
in field theory, see for instance the book of Schwinger 1 1 . 

From this Lagrangian we find the equations of motion and 
the momenta. For instance 



p<p=<p- 



and 



dt 2 



dtp 



dA 
~dt 



(35) 



(36) 



These equations linking the momentum p v to the velocity and 
the equation for the acceleration of the variable tp are remark- 
ably similar to those found in electromagnetism if one inter- 
prets A as the vector potential. 

We shall now use the same trick of considering the con- 
strained Lagrangian dynamics of the field to find the cou- 
pled equations for the field and moving particles. We interpret 
the variable as a rotation velocity and E 2 as a potential 
energy. For notational simplicity we consider unit mass parti- 
cles and a system of units where eo = I& — 1. We find the 
following Lagrangian: 



(37) 



rf 3 r A • (E - curl © + 2_, 1i^( r ~ r t)*i) 



Here the Lagrange multiplier A imposes the kinematic con- 
straint eq. Jl 31 in a manner analogous to the rolling constraint 
for the gears. We find the equations of motion by the usual 
variational calculus: it is useful to note that the curl opera- 
tor is self adjoint with appropriate boundary conditions so that 
j A ■ curl B d 3 r = j B ■ curl A d 3 r 



50 : curl A = 

<5E : E + A = 
dA 



(38) 



5r 



fj + 9i-rr-H S rad A) = 



SA : E — curl + ^ <?i<5(r - r i )r l = 

i 

The variation in Sr^ can be rewritten by using the identity 
grad (v. A) = (v. grad )A + v A curl A and by noting that 
! = (JL+ V .grad) 



q t A - qii i A curl A 
: qi(E + r, A curl 0) 







(39) 



These are the normal equations of electromagnetism if we 
identify with B. The Lagrangian corresponds to the tem- 
poral gauge where the scalar potential (j) = 0. A gauge trans- 
formation A — * A + grad ^(t) generates additional terms in 
the Lagrangian of the form </>( div E — p) with </> = vp. 



We can eliminate from the Lagrangian via the Thomson- 
Routh treatment of kinesthenic variables: Consider the modi- 
fied action C = L — pe®. We find that 



^ 2 



(curl A) 2 i E 2 
2 



C = > ; ^- / d 3 r f v ^~ y + ^- ) + / d 3 r A-(E+J) 



(40) 

which we recognize 1 1 as a more conventional Lagrangian for 
electrodynamic systems. 

We construct the Hamiltonian using Dirac's procedure with 
two constraints: 



(Pi - qiA(vi)) 2 f (p e + curl A) 



MP a + 7(Pe - A) d 6 r 



(41) 



where \i, 7 are the multipliers for the primary constraints. The 
initial conditions are p# = A, p^ = and pg = which are 
conserved by the equations of motion in the same way that 
Gauss's law is conserved in the Monte-Carlo formulation. On 
the physically relevant surface the constraint terms are identi- 
cally zero; the extended Hamiltonian still has the normal in- 
terpretation as the conserved total energy. 

Such a description of the electromagnetic field in terms of 
rotors was known to FitzGerald in the nineteenth century as a 
mechanical analogy 1213 of the ether. A square array of wheels 
was constructed; neighboring wheels were connected by an 
elastic band. When two neighboring wheels turn at the same 
angular velocity the elastic band has constant length and the 
elastic energy is constant. When there is a difference of rota- 
tional velocity between wheels the elastic energy of the bands 
changes. Assuming linear elasticity for the elastic bands one 
finds an exact mapping of electromagnetism onto a mechani- 
cal problem. This model was most important in the history of 
electromagnetism: FitzGerald used this model in the very first 
calculation of radiated power from moving charges. 



C. Statistical Mechanics 

The interpretation of as the angular velocity of a rotor 
suggests that it could be coupled to a thermostat to improve 
equilibration of the field degrees of freedom; the linear equa- 
tions that we have found for the electric field are likely to equi- 
librate rather slowly. If we add coupling to external noise, £ 
and friction, T, in eq. (I38i we find 



= curl A - T0 + C(t) 



(42) 



thus coupling the angular velocity to an arbitrary thermostat 
leads to violation of the Maxwell equation div B = 0. 
The partition function is calculated from 



Vp Pq e 



-0H 



(43) 



The integral is over the canonical coordinates q and momenta 
p. The integration region is the set of configurations available 
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to the equations of motion. We thus implicity include delta 
function constraints on and p^. Integration over the mo- 
menta is easy to perform in the presence of Langevin noise 
which destroys the constraints and conservation laws associ- 
ated with the variable pg in Maxwell's equations. What re- 
mains is the integral over the electric fields and particle posi- 
tions. If the dynamics were ergodic we would integrate over 
all values of the field. However Maxwell's equations, even in 
the presence of noise on the momentum degree of freedom, 
include Gauss's law. This constrains the electric field and the 
partition function is given by 

Z c = J Vn J VE e- / ^ d " r J| S( div E - p/e ) 

(44) 

where, now, all degrees of freedom are freely integrated over. 
It is this constrained configurational integral 7 that leads to ef- 
fective Coulomb interactions. 

Combining eqs. (II 31 and 1421 we find the equation for the 
electric field: 

* ( * + r) E = V^E- grad P - (v + |) J + curl ft) 

(45) 

In the limit of low frequencies we can ignore compared to 
r and find an equation entirely equivalent to eq. Jl 51 . The 
damping strongly modifies the large scale nature of the elec- 
tric field dynamics. 

This result seems quite remarkable. It is known from the 
work of Heaviside 14 that the electric field of a moving par- 
ticle is strongly modified at velocities which approach c, the 



speed of light. Despite this eq. J44i implies that the average 
interaction between particles is independent of this longitudi- 
nal contraction of the electric field. Either this is consequence 
of the full Maxwell equations which has not yet been or ex- 
plored or it is a consequence of relaxing the constraint on the 
divergence of the magnetic field leading to the modified large 
scale properties implied in eq. ( 145 1 . We leave further study of 
this problem to a future publication. In either case this could 
permit study of the Coulomb interacting particles via direct 
integration of the Maxwell equations in molecular dynamics 
simulations. 



V. CONCLUSIONS 

We have analyzed a Monte-Carlo algorithm for the simu- 
lation of long ranged Coulomb interactions. We have seen 
that propagation of the electric degrees of freedom is diffu- 
sive. By construction the dynamics sample the equilibrium 
Boltzmann distribution of the charged system. The locality 
of the algorithm allows fast and simple implementations even 
on multiprocessor computers with high communication over- 
heads. We have verified that the Monte-Carlo algorithm re- 
produces well known features of the two component plasma 
such as screening. Our law for the local update of the electric 
field after movement of a particle, fig (1), is a discretized ver- 
sion of the Maxwell displacement current. The algorithm has 
been shown to be closely related to mechanical models of the 
ether introduced in the generation that followed Maxwell. 
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